{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. Reactor\n",
    "\n",
    "(c) 2019, Dr. Ramil Nugmanov; Dr. Timur Madzhidov; Ravil Mukhametgaleev\n",
    "\n",
    "Installation instructions of CGRtools package information and tutorial's files see on `https://github.com/cimm-kzn/CGRtools`\n",
    "\n",
    "NOTE: Tutorial should be performed sequentially from the start. Random cell running will lead to unexpected results. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pkg_resources\n",
    "if pkg_resources.get_distribution('CGRtools').version.split('.')[:2] != ['3', '1']:\n",
    "    print('WARNING. Tutorial was tested on 3.1 version of CGRtools')\n",
    "else:\n",
    "    print('Welcome!')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data for tutorial\n",
    "from pickle import load\n",
    "from traceback import format_exc\n",
    "\n",
    "with open('reactions.dat', 'rb') as f:\n",
    "    reactions = load(f) # list of ReactionContainer objects\n",
    "\n",
    "r1 = reactions[0] # reaction\n",
    "m6 = r1.reactants[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Reactor* objects stores single transformation and can apply it to molecules or CGRs.\n",
    "\n",
    "Transformations is ReactionContainer object which in reactant side consist of query for matching group and in product side patch for updating matched atoms and bonds "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from CGRtools import CGRreactor   # import of Reactor\n",
    "from CGRtools.containers import * # import of required objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.1. Products generation\n",
    "Reactor works similar to ChemAxon Reactions enumeration.\n",
    "\n",
    "Example here presents application of it to create esters from acids.\n",
    "\n",
    "First we need to construct carboxy group matcher query. Then, ether group need to be specified. \n",
    "\n",
    "Atom numbers in query and patch should be mapped to each other. The same atoms should have same numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acid = QueryContainer()                       # this query matches acids. Use construction possibilities.\n",
    "acid.add_atom({'element': 'C', 'neighbors': 3})   # add carboxyl carbon. Hybridization is irrelevant here\n",
    "acid.add_atom({'element': 'O', 'neighbors': 1})   # add hydroxyl oxygen. Hybridization is irrelevant here \n",
    "acid.add_atom('O')                                # add carbonyl oxygen. Number of neighbors is irrelevant here.\n",
    "acid.add_bond(1, 2, 1) # create single bond between carbon and hydroxyl oxygen\n",
    "acid.add_bond(1, 3, 2) # create double bond\n",
    "print(acid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "methyl_ester = MoleculeContainer()  # create patch - how acrboxyl group should be changed. We write methylated group\n",
    "methyl_ester.add_atom('C', 1) # second argument is predefined atom mapping. Notice that mapping corresponds...  \n",
    "methyl_ester.add_atom('O', 2) # ... to order in already created acid group. Atom 2 is released water.\n",
    "methyl_ester.add_atom('O', 4)\n",
    "methyl_ester.add_atom('O', 3)\n",
    "methyl_ester.add_atom('C', 5)\n",
    "methyl_ester.add_bond(1, 4, 1)\n",
    "methyl_ester.add_bond(1, 3, 2)\n",
    "methyl_ester.add_bond(4, 5, 1)\n",
    "# No bond between atom 1 and atom 2. This bond will be broken. \n",
    "methyl_ester.calculate2d()\n",
    "methyl_ester"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m6.reset_query_marks() # required for correct matching\n",
    "m6 # acid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "template = ReactionContainer([acid], [methyl_ester]) # merge query and patch in template, which is ReactionContainer\n",
    "reactor = CGRreactor(template)                        # Reactor is initialized\n",
    "reacted_acid = reactor(m6)                            # application of Reactor to molecule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reacted_acid.calculate2d(scale=2) # calculate coordinates\n",
    "reacted_acid       # desired methylated ester have been generated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can notice presence of separate oxygen (water) and ester group.\n",
    "\n",
    "The second group can substituted by calling reactor on observed product."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reacted_acid.reset_query_marks() # this is new molecule and query marks need to be set\n",
    "second_stage = reactor(reacted_acid) # apply transformation on product of previous transformation\n",
    "second_stage.calculate2d(scale=2) #  recalculate coordinates for correct drawing\n",
    "second_stage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "second_stage has 3 components in a single MoleculeContainer object. We can split it into individual molecules and place all molecules into ReactionContainer object. Since in CGRtools atom-to-atom mapping corresponds to numbering of atoms in molecules, the resulting product has AAM according to the rule applied. Thus, reaction has correct AAM and nothing special should be made to keep or find it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "products = second_stage.split() # split product into individual molecules\n",
    "react = ReactionContainer([m6], products) # unite reagent and product into reaction. \n",
    "react.fix_positions()\n",
    "react"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For multicomponent reactions one can merge molecules of reactants into single MoleculeContainer object and apply reactor on it.\n",
    "\n",
    "It is possible to generate all available products in case that molecule has several groups matching the query.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m6copy = m6.copy() # let's try to use molecule with several groups mathcing query\n",
    "m6copy.atom(5).isotope = 13 # isotope mark is added to see the difference in products\n",
    "m6copy.reset_query_marks() # query marks need to be calculated\n",
    "enums = set()              # the set enums is used to select structurally diverse products\n",
    "for m in reactor(m6copy, limit=0): # limit=0 is enumeration of all possible products by reactor\n",
    "    print(m)                       # print signatures for observed molecules. Notice presence of water as component of product\n",
    "    m.calculate2d(scale=2)         # recalculate coordinates\n",
    "    enums.update(m.split())        # split product into separate molecules\n",
    "enums = list(enums)                # set of all resulting molecules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m6copy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a look at molecules in set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "enums[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "enums[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "enums[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.2. MetaReactions (reactions on CGRs).\n",
    "Reactor could be applied to CGR to introduce changes into reaction. \n",
    "\n",
    "### 6.2.1. Example of atom-to-atom mapping fixing. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reactions[1] # reaction under study"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cgr = ~reactions[1] # generate reaction CGR\n",
    "print(cgr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cgr.centers_list # reaction has two reaction centers. [10,11,12] - pseudo reaction appeared due to AAM error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "reaction has AAM error in nitro-group\n",
    "\n",
    "![error.png](error.png)\n",
    "\n",
    "Lets try to use Reactor for AAM fixing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nitro = QueryCGRContainer() # construct query for invalid reaction center - CGR of wrongly mapped nitro-group\n",
    "nitro.add_atom({'element': 'N', 'charge': 1, 'p_charge': 1}) # atom 1\n",
    "nitro.add_atom({'element': 'O', 'charge': 0, 'p_charge': -1}) # atom 2. Notice that due to AAM error charge was changed\n",
    "nitro.add_atom({'element': 'O', 'charge': -1, 'p_charge': 0}) # atom 3. Notice that due to AAM error charge was changed\n",
    "nitro.add_atom('C') # atom 4\n",
    "\n",
    "nitro.add_bond(1, 2, {'order': 2, 'p_order': 1}) # bond between atoms 1 and 2. Due to AAM error bond is dynamic ('2>1' type) \n",
    "nitro.add_bond(1, 3, {'order': 1, 'p_order': 2}) # bond between atoms 1 and 3. Due to AAM error bond is dynamic ('1>2' type) \n",
    "nitro.add_bond(1, 4, 1) # ordinary bond\n",
    "print(nitro)\n",
    "# this query matches reaction center in CGR appeared due to AAM error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nitro < cgr # query matches CGR of reaction with error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "valid_nitro = MoleculeContainer() # construct nitro group without dynamic atoms. Notice that atom order should correspond object nitro\n",
    "valid_nitro.add_atom({'element': 'N', 'charge': 1}) # ordinary N atom\n",
    "valid_nitro.add_atom({'element': 'O', 'charge': -1}) # ordinary negatively charged oxygen atom\n",
    "valid_nitro.add_atom('O')                            # ordinary oxygen atom\n",
    "\n",
    "valid_nitro.add_bond(1, 2, 1) # ordinary single bond\n",
    "valid_nitro.add_bond(1, 3, 2) # ordinary double bond\n",
    "print(valid_nitro)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "valid_nitro.calculate2d()\n",
    "valid_nitro # this is correct representation of group."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now time to prepare and apply **Template** to CGR based on reaction with incorrect AAM.\n",
    "\n",
    "Template is Reaction container with query in reactants and patch in products"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "template = ReactionContainer([nitro], [valid_nitro]) # template shows how wrong part of CGR is transformed into correct one.\n",
    "print(template) # notice complex structure of query: CGR signature is given in braces, then >> and molecule signature"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Reactor` class accept single template. Existence of dynamic bond in it is not a problem.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reactor = CGRreactor(template)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Reactor` object is callable and accept as argument molecule or CGR.\n",
    "\n",
    "NOTE: `fixed` is new CGR object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fixed = reactor(cgr) # fix CGR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CGRreactor returns None if template could not be applied, otherwise patched structure is returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(fixed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`C-C(-I).O-C:1:C:C:C(:C:C:1)-[N+](-[O-])=O`  \n",
    "`C-C(.I)-O-C:1:C:C:C(:C:C:1)-[N+](-[O-])=O`\n",
    "\n",
    "One can see that nitro group has no dynamic bonds any more. CGR corresponds only to substitution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fixed.centers_list # reaction center appeared due to AAM error before does not exist. Only 1 reaction center is found"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is depiction of observed CGR (external software was used). Notice absence of wrong reaction center.\n",
    "![cgr3.png](cgr3.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.3.2 Reaction transformation\n",
    "Example of E2 to SN2 transformation.\n",
    "\n",
    "E2 and SN2 are concurrent reactions.\n",
    "We can easily change reaction center of E2 reaction to SN2. It could be achieved by substitution of reaction center corresponding to double bond formation in E2 reaction by the one corresponding to formation of new single bond with base as in SN2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from CGRtools import CGRreactor, CGRpreparer\n",
    "from CGRtools.containers import QueryCGRContainer, ReactionContainer\n",
    "from CGRtools.files import MRVread, SDFwrite\n",
    "from io import StringIO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "e2 = next(MRVread('e2.mrv')) # read E2 reaction from ChemAxon MRV file\n",
    "e2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create CGR query for E2 reaction side\n",
    "e2query = QueryCGRContainer() \n",
    "e2query.add_atom('C', 1) # create carbon with mapping number 1\n",
    "e2query.add_atom('C', 2) # create carbon with mapping number 2\n",
    "# addition of any halogen atom\n",
    "e2query.add_atom({'element': ['I', 'Cl', 'Br'], 'neighbors': 1, 'p_neighbors': 0, 'charge': 0, 'p_charge': -1}, 3)\n",
    "# addition of OH-, RO-, SH-, RS- groups\n",
    "e2query.add_atom({'element': ['O', 'S'], 'neighbors': [0, 1], 'p_neighbors': [0, 1], 'charge': -1, 'p_charge': 0}, 4)\n",
    "\n",
    "e2query.add_bond(1, 2, {'order': 1, 'p_order': 2}) # bond between two carbon corresponds to formation of double from single\n",
    "e2query.add_bond(1, 3, {'order': 1, 'p_order': None}) # bond between carbon and halogen breaks in E2 reaction\n",
    "print(e2query) # it is CGR of E2 reaction center"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "e2_cgr = ~ e2 # compose reaction into CGR\n",
    "e2_cgr.reset_query_marks() # prepare to isomorphism"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CGR is the following (depicted by external software)\n",
    "![cgr4.png](cgr4.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "e2query < e2_cgr # E2 CGR pattern works!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create patch creating SN2 reaction. Notice that ordering of atoms correspond to that of E2 CGR query\n",
    "sn2patch = QueryCGRContainer()\n",
    "sn2patch.add_atom({}, 1) # save atom unchanged. We don't specify atom type since it will be taken from E2 query\n",
    "sn2patch.add_atom({}, 2) # it is central atom. We don't specify atom type since it will be taken from E2 query\n",
    "sn2patch.add_atom({'charge': 0, 'p_charge': -1}, 3) # elements list with same order of elements [I, Cl, Br] could be used as well\n",
    "sn2patch.add_atom({'charge': -1, 'p_charge': 0}, 4)\n",
    "\n",
    "sn2patch.add_bond(1, 2, {'order': 1, 'p_order': 1}) # set carbon - carbon single bond that is unchanged in SN2 reaction\n",
    "sn2patch.add_bond(1, 3, {'order': 1, 'p_order': None}) # this bond is broken in SN2 reaction\n",
    "sn2patch.add_bond(1, 4, {'order': None, 'p_order': 1}) # it corresponds to formation of bond O(S)-C bond in SN2 reaction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reactor = CGRreactor(ReactionContainer([e2query], [sn2patch])) # create template and pass it to Reactor\n",
    "sn2_cgr = reactor(e2_cgr) # apply Reactor on E2 reaction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sn2_cgr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is depiction of CGR produced by Reactor (external software is used)\n",
    "![cgr5.png](cgr5.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# decompose CGR into reaction\n",
    "preparer = CGRpreparer() \n",
    "sn2 = preparer.decompose(sn2_cgr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sn2.calculate2d()\n",
    "sn2 # reaction has the same reagents like E2 above, but products correspond to SN2 reaction"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "cgrtools",
   "language": "python",
   "name": "cgrtools"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
